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A semiclassical (SC) approach is developed for nonequilibrium quantum transport in molecular junctions. 
Following the early work of Miller and White [J. Chem. Phys. 84, 5059 (1986)], the many-electron Hamil- 
tonian in second quantization is mapped onto a classical model that preserves the fermionic character of 
electrons. The resulting classical electronic Hamiltonian allows for real-time molecular dynamics simulations 
of the many-body problem from an uncorrelated initial state to the steady state. Comparisons with exact 
results generated for the resonant level model reveal that a semiclassical treatment of transport provides a 
quantitative description of the dynamics at all relevant timescales for a wide range of bias and gate potentials, 
and for different temperatures. The approach opens a door to treating nontrivial quantum transport problems 
that remain far from the reach of fully quantum methodologies. 



I. INTRODUCTION 

Molecular electronics^ has provided means to study 
the dynamics of open quantum systems, in which one 
considers a small, strongly interacting and highly corre- 
lated region (the molecule and its closest vicinity) cou- 
pled to several large, noninteracting baths (representing 
the fermionic leads and environment). While the equilib- 
rium nature of quantum dynamics in condensed phases 
has been mostly resolved, 21 the intrinsic nonequilibrium 
nature of transport through molecular junctions, along 
with the necessity to treat fermionic degrees of freedom, 
poses a much greater theoretical challenge, and thus re- 
mains poorly understood. 

Several different paths have been taken to improve 
the standard Landauer-Biittiker approach 3 ^ and its gen- 
eralization to the multichannel caseP in order to ac- 
count for electron-electron and electron-phonon correla- 
tions in molecular junctions. The different approaches 
can be classified as perturbative treatments, among 
which the most notable examples use the nonequilibrium 
Green's function (NEGF) formalism]^ and numerically 
exact techniques, perhaps the most prominent of which 
have been time-dependent numerical renormalization 
group techniques'^ anc j th e promising diagrammatic 
approaches based on path integral formulations CUED 

These approaches have been applied to a variety of 
physically interesting problems including the descrip- 
tion of Coulomb and Franck-Condon blockadep^H^ the 
nonequilibrium Kon do problem] 20 ! 21 ! and inelastic elec- 
tron tunneling ! 11 ! 22 ! While successful to a large extent, 
these approaches suffer from several limitations, includ- 
ing difficulties that arise in the treatment of more com- 
plex environments or when inelastic scattering is gov- 
erned by interactions of electrons with soft modes that 



are dominated by large anharmonicities. 

Parallel to these developments, a completely different 
paradigm has been devised based on semiclassical ap- 
proaches. These have provided a useful tool to simu- 
late the dynamics of molecular subsystems coupled to a 
fluctuating environment with significant anharmonicities. 
The most appealing semiclassical treatments have relied 
on the mixed quantum-classical approach^!!!!] anc j on t ne 
semiclassical initial value representation!^!!!] The latter 
have been applied to a variety of physically interesting 
condensed phase problems with remarkable successPSHI] 
It is interesting to note that such approaches have not 
received any attention in the context of nonequilibrium 
quantum t ranspor t , despite being exact in the harmonic 
boson case! 27 ' 42 ' 43 ! 

A major goal of the present work is to show how a 
semiclassical (SC) approach can treat the dynamics of 
a many-body quantum system driven away from equi- 
librium by the application of a bias voltage. We de- 
scribe the transport problem in second quantization, par- 
titioning the space into an interacting region describ- 
ing the molecule and its closest vicinity and a nonin- 
teracting region representing the leads and the environ- 
ment. The approach is based on an SC model for the 
general second-quantized many-electron Hamiltonian by 
Miller and White (MW) , 44 which followed earlier work of 
McCurdy, Meyer, and Miller (MMM)PEH1 on construct- 
ing classical models for electronic degrees of freedom. 
The essence of MWs model is that each fermionic de- 
gree of freedom (i.e., each pair of one-particle annihila- 
tion/creation operators) is described by a classical de- 
gree of freedom (pair of action-angle variables) , while in 
MMM's earlier work each electronic state is described by 
a classical degree of freedom. MWs model is thus a more 
"microscopic" description of the electronic degrees of free- 
dom, and more importantly much more "efficient" (i.e., 
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involving many fewer classical degrees of freedom) for 
many electron systems (where the number of electronic 
states can be much larger than the number of one-particle 
annihilation / creation operators) . 

To assess the accuracy of the proposed approach, we 
focus on the resonant level model and derive working 
expressions to simulate the left, right and total current 
using classical trajectories with quasi-classical initial con- 
ditions. Remarkably excellent agreement in comparison 
to exact results is achieved for a wide range of bias and 
gate voltages and for different temperatures. Our ap- 
proach provides a natural framework to study more com- 
plex molecular transport problems with promising perfor- 
mance. 

The paper is organized as follows. In Section [IT] we 
summarize the semiclassical procedure for constructing a 
classical model of a second quantized Hamiltonian. Sec- 
tion [TllJ describes the resonant level model, which is used 
as a test case to assess the accuracy of the semiclassical 
approach. Exact quantum mechanical results for the res- 
onant level model are provided in Section IV Section [V] 
summarizes the main results and provides a detailed com- 
parison between the semiclassical approach and the exact 
quantum mechanical treatment. The comparisons cover 
a wide range of gate and bias potentials from high to low 
temperatures. Section [VI| summarizes and concludes. 



II. SEMICLASSICAL APPROACH 



A. Mapping 



The approach that has been usecffiEiHH] t construct 
semiclassical models for electronic degrees of freedom is 
to invert (as meaningfully as possible) the Heisenberg 
correspondence relation, 



A 



(2tt 



,-F 



2 k 



e-^ n '- n ^A cl (n,q)dq (1) 



where n = g(n' + n), F is the number of degrees of free- 
dom, and A c \(n,q) is a function of the classical action- 
angle variables. The Heisenberg correpondence relation 
was originally used to obtain approximate matrix ele- 
ments for A from the corresponding classical function of 
action-angle variables; the semiclassical goal here is to 
obtain a classical function of action-angle variables that 
corresponds (as best as possible) to the given quantum 
mechanical matrix elements of operator A. The formal 
inverse of the Fourier transform in Eq. (JTJ gives the angle 
dependence of the classical function as a Fourier series, 



A c \{n,q) = ^2 
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Here, however, the procedure becomes nonunique (and 
nonexact) because the classical action variable n has to 
interpolate for all real numbers, while the matrix ele- 
ments are only defined for integer values of n ± | . 



One way to proceed is to use the spin-matrixmapping 
(SMM) method of Meyer and Miller (MM),™ 3 wn i c h 
utilizes the fact that a general two-state system is equiv- 
alent to a spin | system. Thus any 2x2 matrix A n , n r, 
7i, n' G {0, 1}, can be written as a linear combination of 
the three spin matrices S x , S y , and S z (and the 2x2 
identity matrix). One then uses the classical expressions 
for the spin angular momentum, 



S x — \J a 2 — m 2 cos(g) 

S y = \] <j 2 — m 2 sin(g) 
S z = m. 



(3a) 

(3b) 
(3c) 



in terms of the action-angle variables (m, q). The quan- 
tum values of the action variable m (the projection quan- 
tum number) are ±~, but it is convenient to have the 
quantum values of the action variables be (unoccupied) 
and 1 (occupied). The elementary canonical transforma- 
tion ri = m + h makes this change, giving the classical 
function of action-angles (n, q) for a general 2x2 matrix 
A as 



A(n,q) = (l-n)A Q0 + nA 
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+ - n 2 + A (A 10 e l * + A 01 e^ 9 ) (4) 

where A = a 2 — 4. The value of a 2 will be discussed 
below. 

To employ this approach on a given second-quantized 
Hamiltonian, MW used this SMM model for each term in 
a Hamiltonian. The Hamiltonian in the present work has 



terms of the form a\ hi and a] &j . The matrix elements of 
the operator a\a,i are diagonal, 



a]a.i 



J[Sn 



(5a) 



which corresponds to Eq. Q with A ao = and An = 1 
for the ith degree of freedom, and to the identity matrix 
for all other degrees of freedom. The SMM model thus 



gives 



a ' a,- i — y = 



The matrix elements for a\a,j are 



(5b) 
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where we have assumed that i < j. In the case i > j, the 
product over p goes from j + 1 to i — 1: it always includes 
the states between (but not including) i and j. That 
product is the result of the anticommutation relation of 
fermionic creation/annihilation operators. 

These matrix elements are separable products of 2 x 2 
matrices: for the ith degree of freedom (corresponding to 
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the factor 5 n > yn '+i) the 2x2 matrix has Aio — 1 as the 
only nonzero matrix element, so that Eq. Q gives the 
corresponding classical function as 



- nj + A e 
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(6b) 



For the jth degree of freedom (i.e., the factor £ n '., n ._i), 
the only nonzero matrix element is Ao± = 1, resulting in 
the classical function 



nj + A e 



"ij 



(6c) 



For each degree of freedom p. with p between i and j 
(in normal ordered form), the factor S n ^ np (— l) ,lp from 
Eq. ( 6a I corresponds to a diagonal 2x2 matrix in Eq. Q 
with% = 1 and An = -1. This SMM method thus 
gives the classical function 



1 — 2n v 



(6d) 



The net result is 



+ A) ( nj — n 2 + \) e l ^~ q ^ 



Y[fb(n p ) 



(7) 



where the SMM method gives fb(n p ) = 1 — 2n p . For 
the present application, however, we found that using 
these factors /b(np) significantly underestimates the cur- 
rent (see section [VJ, while excellent results were obtained 
by omitting these factors, i.e., by setting fb{n p ) = 1. We 
also tried an alternative, fb{n p ) — exp(i7m p ), but it had 
the same defect. The issue seems to be that with the 
quantum values of n p 6 {0, 1}, the product of these fac- 
tors is ±1, while the classical value of n p frequently leads 
to fb(n p ) < 1. The product of many such factors leads to 
factors much less than 1, and thus off-diagonal coupling 
that are much too small. 

We now return to the value of a (which determines 
A). For a classical spin with s — |, a 2 = s 2 = j, such 
that A = 0. However, MM suggested using the Langer 
modified value a 2 = (s + |) = 1, giving A = |. Here, 
we suggest instead using the quantum value of a 2 = s(s + 
1) = |, and thus A = \. The results in section |v| give 
empirical reasons for this choice. 

The choice A = \ also can be justified with the follow- 
ing argument. Because the goal is to generate a mapping 
that accounts for the correct dynamics, it is reasonable to 
try to match the short-time dynamics of a simple prob- 
lem. Consider the Hamiltonian H = k ^a^a_B + flgaij , 
with initial population in state A. The first time deriva- 
tive of the initial (t = 0) quantum population of state 
A is zero, as is that of our semiclassical model (aver- 
aged over angles). The second derivative of the initial 
quantum population of state A is —2k 2 , while the semi- 
classical model gives — 4fc 2 A. If these are to be equal, 
then we must set A = 



To summarize, the modifications we have made to the 
original MW model are (1) to use a value of A = | rather 
than the previously suggested Langer modified value of 
|, and (2) to set the factors fb(n p ) in equation ^ to 
unity. 



B. Intial conditions 

The semiclassical mapping also requires the selection 
of initial conditions for the action-angle variables. It 
is clear that a naive approach based on the classical 
mapping described above will not provide the correct 
statistical treatment: Even for the one-particle system 
H = a^a i— > n, it does not reproduce the quantum parti- 
tion function when the action-angle variables are sampled 
from the corresponding classical thermal distribution. 

To recover the correct statistical behavior (at least at 
time t = 0) we use a quasi-classical procedure. Since we 
are interested in a noncorrelated initial state with ther- 
mally populated leads and an unpopulated quantum dot, 
we can populate each degree of freedom independently. 
We enforce quantum statistics on the initial conditions 
for each degree of freedom i by setting the initial action 
rii to either or 1 such that the expectation value of the 
action (rii), averaged over the set of initial conditions, 
satisfies the Fermi distribution, 



/(e* - Mi) 



(8) 



where fii is the chemical potential of the lead in which 
mode i is located and j3 = 1/T is the inverse temperature. 
The angle variable q t is selected at random between and 
2tt. 



III. MODEL HAMILTONIAN 

We use the resonant level (Landauer) model as an ex- 
ample of quantum transport. It consists of a single quan- 
tum dot state coupled to two electrodes (left and right) 
according to the Hamiltonian 



H = e aja 



N 



N 



k=l 



k=l 



a\a 



(9) 



where eo is the energy of the isolated quantum dot (and 
will also be used to model a gate voltage), is the energy 
associated with the electrode mode fc, and is the cou- 
pling between the quantum dot and the electrode mode 
k. 

Using the procedure described in section [H] we con- 
struct the classical model of this Hamiltonian in action- 
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angle variables: 



N N I ~ 

H(n, q) = e n + 2J t k n k +^Jk \ n - n\ + - 



fc=i 



x \ n k - ni 



fc=i 



J(qo-qk) 



N 



Y[f b (n p )+R, 



P =i 



(10) 



The left current is given by the change in occupancy of 
the left electrode, with the right current defined analo- 
gously and the total current given by half the difference of 
the two. For the Hamiltonian in Eq. ^ , the Heisenberg 
time derivative gives the left current as: 



d 

dT 



< — { ^a\a k 

\k£L I 



(11) 
(12) 



\keL 



where L is the set of states in the left electrode. 

We can choose to take the semiclassical approximation 
before the time derivative by mapping the occupation 
to a classical quantity and taking its time derivative, or 
we can take the result of the quantum time derivative 
and map that to a classical quantity. For this system, 
our mapping method gives formally equivalent results for 
either procedure. 

The electrodes are described within the wide band 
limit with a sharp cutoff at high and low energy values: 



Jl/r{z) = 



L/R 



(13) 



(l + e^-f)) (l + e-^+f)) 
where, in all results reported below we use T l = Y R = g 



r = T L +T R , A = 5r, and B = 20T. For the semiclassical 
mapping, we use a uniform discretization to select the 
energies of the leads' states, and thus the couplings are 
given by 



ife(efe) = 



J(e fc )Ae 
2tt 



(14) 



IV. EXACT QUANTUM MECHANICS 

An exact quantum mechanical solution for the tran- 
sient current for the Hamiltonian specified above is 
straightforward to derive, and thus provides means to as- 
sess the accuracy of the semiclassical treatment. Under 
the assumptions of no correlation at t = 



4 (0) ao (0)\ = (4 (OK (0)\ =0 



(15) 



and a Boltzmann distribution for the leads' populations 

(4 (0) a k , (0)) = / (e fc - t* LtR ) 5 kk >, (16) 
where / (e) is the Fermi function and as before hl,r is 
the chemical potential for the left (L) or right (R) lead, 
one can derive an exact expression for the left current 



given in Eq. ( 12 1 



l J L (w')dw' 



for any given initial dot population 

'4(0)a (0)\=n G [0,1]. 



(17) 



(18) 



In the above, Jl (u/) is given by: 



Jl («;') 



£ p i i (-^<n + (-|s<(^-^)) f ) 

2tt7 — I -iu)' (w - f - ±X L , R (w)) 
(S £ (u-u')^n (0) 



in 
n 



1 (£ £ ( W - c')) f (-f E< R ( W ) + (-|S< R ( W - a/)) 

- ^ - f - H (^,« (« - ^')) f ) - f - (a;)) 



dcj, 



(19) 



£ (cj) is the self energy, and £< (w) is the lesser self en- 
ergy, both specified below. A similar expression can be 
derived for the right current (I R (t)) by the replacement 
L <H- R. The total current is given by the different of the 



left and right currents, I (t) = fL(t) ' fa( ' ) . 

As a check on the above, we discuss two known limits 
for the current. In the limit t — > 0, it is simple to show 
that the current vanishes. This is a result of the initial 
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preparation of an uncorrelated state and can be derived 
with the help of the initial value theorem 

lim I L (t) = ^Im { lim -iuj 1 J {ui') 1=0, (20) 

where we assumed that the self-energies vanish at the 
boundaries in the frequency domain: 



To recover the Landauer expression for the current at 
steady state^^we take the limit t —> oo in Eq. (17 1, this 
time with the aid of the final value theorem 



lim E (oj) — lim E (oj) 



0. 



(21) 



2e 

lim (t) — ^Im 



lim —ioj'J(oj') 

-iuj'— f 



1 



h 2tt 



x (s L H) t (|E^M + ( 



do; 



(22) 



Rearranging Eq. (22 I and using the well known represen- V. RESULTS 



tation of the self energies in terms of the real (^l,r (oj)) 
and imaginary (Tl^r (oj)) portions: 

E L , R (oj) = A L , R (oj) - \%Tl,r (w) (23) 



s l r ( w ) = (J-l,r) Tl,_r M 



(24) 



we finally arrive at the Landauer expression for the cur- 
rent: 

/>00 



(f(hw - n L ) - f(hw - n R )) 



Tr (oj) T l (oj) 



-do: 



where the rescaled energy e~ is given by: 
e~o = £o + Al M + A fl (w) 

andr(w) = r L (w) + r fl (w). 

We choose 



? l /r(u) = 



L/R 



(25) 



(26) 



(27) 



where Tl = Tr = h. A and B are defined in Eq. (13 1 
K l /r(oj) is obtained from the Kramers-Kronig relation 



Al/,rH 



1 



-PP 



da; 



7T ./_~, 0/ — OJ 

where pp denotes the Cauchy principle value. 



(28) 



Our simulations used the semiclassical Hamiltonian of 
Eq. ( |10[ ) in action-angle variables, with initial conditions 
selected as described in section |TT] Numerical integra- 
tion of the trajectories was performed with the sixth- 
order Gear predictor-corrector algorithm, modified to 
allow adaptive timesteps (necessary due to the square 
roots in the Hamiltonian). The maximum step size was 
At = 0.01/i/r. We used 400 states per lead and reported 
results with 2 x 10 5 trajectories (see below for further dis- 
cussion on these points). A calculation to time t = 8h/T 
with 400 modes per electrode and 2 x 10 5 trajectories 
required less than one hour walltime with 120 computa- 
tional cores. 

In Figs. 1-3, we plot the left (lower panels), right (mid- 
dle panels) and total (upper panels) currents as a func- 
tion of time for different source-drain voltages (Fig. [TJ, 
gate voltages (Fig. [2]) and temperatures (Fig. [3). The left, 
right and total currents show a distinct time-dependence 
and decay to the same value at steady state. In these 



plots, the exact quantum results of section IV are pre- 
sented as lines, and the results of the semiclassical model 
presented in section |TT] are presented as symbols. 

Figure [T] shows the time-dependent current at a wide 
range of source-drain voltages eV = (J-l — M-Rj with eo = 
(zero gate voltage) and temperature T = T/3. The 
largest values of this bias are limited by the width of our 
band, B = 20F [see Eq. (fl3|]. The most striking result 
shown in Fig. [T| is the excellent agreement between the 
semiclassical method and the exact quantum mechanical 
result. In some cases it is difficult to distinguish the two. 
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Figure 1. Plots of the left (lower panel), right (middle panel) 
and total (upper panel) transient current for different values 
of the source-drain voltage eV = fiL — l-iR- Other model 
parameters are: T = T/3, e = and N L = N R = 400. Solid 
line and symbols correspond to exact quantum mechanical 
and semiclassical results, respectively. 



Figure 2. Plots of the left (lower panel), right (middle panel) 
and total (upper panel) transient current for different values 
of the gate voltage eV g = eo- Other model parameters are: 
Hl = -flR = T, T = T/3 and N L = N R = 400. Solid line 
and symbols correspond to exact quantum mechanical and 
semiclassical results, respectively. 



The semiclassical approach captures both the evolution 
of the current at early and intermediate times, and its 
decay to the correct steady-state value at longer time. 
The agreement is quantitative for all source-drain volt- 
ages studied. The pronounced oscillations in the left and 
right current as they decay to steady state result from 
the finite band width of the leads, and are also captured 
by the semiclassical method. 

At zero source-drain bias, the left and right currents 
show a significant transient effect until the zero-current 
steady-state value is reached, while the total current is 
identically zero at all times. This suggests that for other 
systems the total current would be a better observable 
since one can infer the steady-state result from relatively 
shorter timesP^H2l For the largest bias presented, / = 
0.241T, which approaches the infinite bias limit of I — 
0.25r. 

In Fig. [2] we plot the time-dependent current for differ- 
ent gate voltages, eVc- The gate voltage is determined 
by the parameter eVc — eo in the Hamiltonian The 
line for eVa = — T is not plotted for the total current be- 
cause the quantum result overlaps the eVa — T line (the 
semiclassical results for these two series nearly overlap as 
well). For nonzero gate voltages, we find that the semi- 



classical method is not as accurate as for the case where 
eVa = 0, although it still captures the correct trends at 
all times. In particular, the even more pronounced oscil- 
lations and the decay of the amplitude of the oscillations 
are captured by the semiclassical treatment. Addition- 
ally, the steady-state currents for gate voltages of equal 
magnitude but opposite signs (e.g., bVg — ±r) are equal 
(within the numerical noise), just as in the quantum me- 
chanical results for these parameters. 

The effect of temperature is shown in Fig. [3] We ex- 
plored a range temperatures from T/5 to 2r, which cov- 
ers the classical to quantum regimes. We find excellent 
agreement for all temperatures, even when considering 
the left and right currents which exhibit transient phe- 
nomena on longer timescales and are thus more difficult 
to described A closer examination of the results reveals 
that the semiclassical treatment is a bit more accurate 
for higher temperatures (as might be expected), yet the 
overall agreement for all temperatures is quite surpris- 
ing, in particular in view of the failures of semiclassical 
treatments at low temperatures for other systems.^ 

Fig. [3] deals with the choices we made in the semiclas- 
sical procedure described in section [IT] The upper left 
panel displays the effect of choosing different values of 
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Figure 3. Plots of the left (lower panel), right (middle panel) 
and total (upper panel) transient current for different tem- 
peratures. Other model parameters are: = —/-ir = T, 
eo = and Nl = Nr — 400. Solid line and symbols corre- 
spond to exact quantum mechanical and semiclassical results, 
respectively. 



A (or, equivalently, a 2 ). The choice A = h corresponds 
to the quantum a 2 = s(s + 1), and is what was used in 
the results shown in Figs. 1-3. The choice A = I is the 
the Langer-modified a 2 = (s + \) 2 suggested by MW. 
Stock has treated the Langer modification as a free pa- 
rameter,^ and so in that spirit we include the choices 
A = 1 and A — |, the latter being near Stock's result of 
0.3125 (for a different system). In particular, note the 
short-time behavior: The slope of the semiclassical cur- 
rent at time t = can be shown to be proportional to A, 
and A = | best matches the exact quantum result. This 
is consistent with the short-time analysis in section |IT| 

The upper right panel of Fig. [4] shows how the function 
fb{n), which relates to the anticommutation relations, 
affects the results. The SMM form for these factors, 
fb(n) = 1 — 2n, is labeled as "Linear," while the alter- 
nate we explored, fb(n) = exp(i7m), is labeled "Exp." 
The method we used in the results presented in Figs. 1-3 
was to omit these terms (i.e., fb(n) — 1) and is labelled 
"None." Clearly, the linear and exponential versions un- 
derestimate the current (the linear version tending to 
zero steady-state current), while omitting these factors 
gives extremely good agreement. This can be traced to 
the fact that quantum mechanically, the values of these 




Figure 4. Plots of the total transient current for different val- 
ues of the Langer modification (upper left panel), for different 
choices of the factor fb(n) (upper right panel), for different 
numbers of states in the leads (lower left panel), and with 
different numbers of trajectories (lower right panel). In all 
panels hl = — fJ-R = T, e = 0, T = F/3 and Nl = Nr. 



factors are either — 1 or 1 which can be absorbed into the 
coupling constants tk as a phase factor. One can prove 
that the current is independent of the sign of tk for the 
initially uncorrelated case. However, in the semiclassi- 
cal mapping the values of these factors are continuous. 
They can therefore significantly change the effective value 
of the coupling constants, decreasing the current. It is 
also notable that the oscillations in the linear and expo- 
nential versions are more pronounced than when these 
factors are omitted. 

The lower left panel shows the importance of using 
a sufficient number of modes per electrode. Since the 
semiclassical simulations are performed for a finite num- 
ber of states in the electrodes, the "steady-state" current 
we obtain cannot hold to infinite time. All other re- 
sults presented involved 400 modes per electrode, which 
Fig. [4] suggests holds the steady state until at least time 
t ps 5h/T. 

The lower right panel of Fig. [4] shows how the results 
converge with the number of trajectories. Even at 2.4 x 
10 4 trajectories, the approximate result is clear. The 
results presented in Figs. 1-3 are with 2x 10 5 trajectories. 



VI. CONCLUDING REMARKS 

The present classical model for electronic degrees of 
freedom (dofs) is seen to provide an excellent descrip- 
tion of the dynamics of electron transmission in a sim- 
ple model of a molecular transistor driven out of equilib- 
rium. The near quantitative results, obtained over a wide 
range of system parameters (including quite low tempera- 
tures), are especially remarkable since the model is imple- 
mented with our rather crude quasi-classical approach: 
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Quantized initial conditions are selected for the action- 
angle variables, but then the dynamics are completely 
classical (i.e., generated by Hamilton's equations, with- 
out any more sophisticated semiclassical input). One 
is accustomed to quasi-classical treatments being semi- 
quantitative, e.g., for quantized vibrational and rota- 
tional dofs — and even for electronically nonadiabatic 
transitions between two or three electronic states — but 
the level of accuracy seen in the present work is much bet- 
ter than that reported in such cases. It may be that the 
dense set of electronic states that constitute the left and 
right electrodes in the present molecular model leads to 
a rapid decay of the initial factorized state to the steady 
state, and thus makes the overall system behave more 
classically. 

The present results encourage one to expand the treat- 
ment to more realistic molecular models of such phe- 
nomena. It is obvious how to include nuclear degrees of 
freedom in the classical molecular dynamics simulations 
along with the classical model for the electronics dofs. 
Also, MW's classical electronic model is able to include 
electron correlation, i.e., two electron interactions in the 
second-quantized Hamiltonian of the type 



E 

i,j,k,l 



\kl) a\a]aia k , 



(29) 



which would be necessary to describe effects such as 
Coulomb blockade. Though the form of this more gen- 
eral classical electronic model has been given, it has never 
been applied to any nontrivial examples. It will be of in- 
terest to see how well the approach will work for such 
more realistic molecular models. 
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